{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d73019f7",
   "metadata": {},
   "source": [
    "### Inversion example \n",
    "\n",
    "* Variant ID: 605947\n",
    "* ENSG00000114547 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "30669239",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import pysam \n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2dfb0e40",
   "metadata": {},
   "outputs": [],
   "source": [
    "wkdir = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3\"\n",
    "wkdir_path = Path(wkdir)\n",
    "\n",
    "# gene and variant \n",
    "misexp_vrnt_id = \"605947\"\n",
    "misexp_gene_id = \"ENSG00000114547\"\n",
    "out_dir= wkdir_path.joinpath(f\"6_misexp_dissect/other_mechanisms/{misexp_vrnt_id}_{misexp_gene_id}\")\n",
    "out_dir.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "vcf_path=\"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/lof_missense/data/sv_vcf/filtered_merged_gs_svp_10728.vcf.gz\"\n",
    "sv_info_path = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/lof_missense/data/sv_vcf/info_table/final_sites_critical_info_allele.txt\"\n",
    "# inputs\n",
    "star_pass_2_out = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq_n5188_v97/results/star_pass2_2ndpass/\"\n",
    "ge_matrix_flat_path = wkdir_path.joinpath(\"2_misexp_qc/misexp_gene_cov_corr/tpm_zscore_4568smpls_8610genes_flat_misexp_corr_qc.csv\")\n",
    "wgs_rna_paired_smpls_path=wkdir_path.joinpath(\"1_rna_seq_qc/wgs_rna_match/paired_wgs_rna_postqc_prioritise_wgs.tsv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dc3fdf52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load SV info \n",
    "sv_info_df = pd.read_csv(sv_info_path, sep=\"\\t\", dtype={\"plinkID\":str})\n",
    "# get variant ID information \n",
    "vrnt_id_info_df = sv_info_df[sv_info_df.plinkID == misexp_vrnt_id]\n",
    "chrom, start, end = [vrnt_id_info_df[col].item() for col in [\"chr\", \"pos\", \"end\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c61eab0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ge_matrix_flat_df = pd.read_csv(ge_matrix_flat_path, sep=\",\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bafc732d",
   "metadata": {},
   "source": [
    "### Expression of SV carriers and non-carriers "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "93caf1a4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gene IDs passing filters: 8610\n",
      "RNA-seq sample IDs passing QC: 4568\n",
      "Loading input VCF ...\n",
      "VCF loaded.\n",
      "Subset VCF to samples with RNA-seq ...\n",
      "Number of samples in VCF with RNA ID and passing QC: 2640\n",
      "Number of RNA IDs passing QC: 2640\n"
     ]
    }
   ],
   "source": [
    "misexp_exp_df = ge_matrix_flat_df[ge_matrix_flat_df.gene_id == misexp_gene_id]\n",
    "# subset to SV samples \n",
    "gene_id_pass_qc_set = set(ge_matrix_flat_df.gene_id.unique())\n",
    "print(f\"Gene IDs passing filters: {len(gene_id_pass_qc_set)}\")\n",
    "smpl_id_pass_qc_set = set(ge_matrix_flat_df.rna_id.unique())\n",
    "print(f\"RNA-seq sample IDs passing QC: {len(smpl_id_pass_qc_set)}\")\n",
    "# egan ID, RNA ID sample links \n",
    "wgs_rna_paired_smpls_df = pd.read_csv(wgs_rna_paired_smpls_path, sep=\"\\t\")\n",
    "egan_ids_with_rna = wgs_rna_paired_smpls_df[wgs_rna_paired_smpls_df.rna_id.isin(smpl_id_pass_qc_set)].egan_id.tolist()\n",
    "# load VCF and subset to EGAN IDs with RNA \n",
    "print(\"Loading input VCF ...\")\n",
    "vcf_path = vcf_path\n",
    "vcf = pysam.VariantFile(vcf_path, mode = \"r\")\n",
    "print(\"VCF loaded.\")\n",
    "print(\"Subset VCF to samples with RNA-seq ...\")\n",
    "vcf_samples = [sample for sample in vcf.header.samples]\n",
    "vcf_egan_ids_with_rna = set(egan_ids_with_rna).intersection(set(vcf_samples))\n",
    "vcf.subset_samples(vcf_egan_ids_with_rna)\n",
    "vcf_samples_with_rna = [sample for sample in vcf.header.samples]\n",
    "print(f\"Number of samples in VCF with RNA ID and passing QC: {len(vcf_samples_with_rna)}\")\n",
    "# subset egan ID and RNA ID links to samples with SV calls and passing QC \n",
    "wgs_rna_paired_smpls_with_sv_calls_df = wgs_rna_paired_smpls_df[wgs_rna_paired_smpls_df.egan_id.isin(vcf_samples_with_rna)]\n",
    "# write EGAN-RNA ID pairs to file\n",
    "rna_id_pass_qc_sv_calls = wgs_rna_paired_smpls_with_sv_calls_df.rna_id.unique().tolist()\n",
    "print(f\"Number of RNA IDs passing QC: {len(rna_id_pass_qc_sv_calls)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4f9c8b7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_rna_id_carriers(vrnt_id, chrom, start, end, wgs_rna_paired_smpls_df):\n",
    "    chrom_num = chrom.split(\"chr\")[1]\n",
    "    vrnt_gt_egan_dict, count = {}, 0\n",
    "    found_vrnt_id = False\n",
    "    for rec in vcf.fetch(chrom_num, start-1, end): \n",
    "        vcf_vrnt_id = str(rec.id)\n",
    "        if vrnt_id == vcf_vrnt_id:\n",
    "            found_vrnt_id = True \n",
    "            # collect genotypes\n",
    "            gts = [s[\"GT\"] for s in rec.samples.values()]\n",
    "            for i, gt in enumerate(gts): \n",
    "                vrnt_gt_egan_dict[count] = [vrnt_id, vcf_samples_with_rna[i], gt]\n",
    "                count += 1 \n",
    "    if not found_vrnt_id: \n",
    "        raise ValueError(f\"Did not find {vrnt_id} in {vcf_path}\")\n",
    "    vrnt_gt_egan_df = pd.DataFrame.from_dict(vrnt_gt_egan_dict, orient=\"index\", columns=[\"vrnt_id\", \"egan_id\", \"genotype\"])\n",
    "    egan_id_carriers = vrnt_gt_egan_df[vrnt_gt_egan_df.genotype != (0, 0)].egan_id.unique()\n",
    "    rna_id_carriers = wgs_rna_paired_smpls_df[wgs_rna_paired_smpls_df.egan_id.isin(egan_id_carriers)].rna_id.unique()\n",
    "    if len(egan_id_carriers) != len(rna_id_carriers): \n",
    "        raise ValueError(\"Number of RNA and DNA samples does not match.\")\n",
    "    return rna_id_carriers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "09d133cc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of carriers with RNA: 10\n"
     ]
    }
   ],
   "source": [
    "rna_id_carriers = get_rna_id_carriers(misexp_vrnt_id, chrom, start, end, wgs_rna_paired_smpls_df)\n",
    "print(f\"Number of carriers with RNA: {len(rna_id_carriers)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ee8e85b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "misexp_exp_sv_carriers_df = misexp_exp_df[misexp_exp_df.rna_id.isin(rna_id_pass_qc_sv_calls)].copy()\n",
    "misexp_exp_sv_carriers_df[\"status\"] = np.where(misexp_exp_sv_carriers_df.rna_id.isin(rna_id_carriers), \"Carrier\", \"Non-carrier\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "14acb1e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "misexp_exp_sv_carriers_df[\"misexp\"] = np.where((misexp_exp_sv_carriers_df.TPM > 0.5) & (misexp_exp_sv_carriers_df[\"z-score\"] > 2), \"Misexpressed\", \"Control\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "86c05fd5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAAHqCAYAAAC9YcJAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeVxV1f7/8fcBmRTBAcUJp5xQnL2ZmEOpqHlNy5KGr3MDappzeb051e/aZJmlVKbNqZXDrSSHnGcTJU3IWTHDARVwAgT27w+upwjQc/AcDltfz8fjPGKvvfben0M98u3ae61tMQzDEAAAAIo8N1cXAAAAANsQ3AAAAEyC4AYAAGASBDcAAACTILgBAACYBMENAADAJAhuAAAAJkFwAwAAMIk7LrgZhqGUlBSx7jAAADCbOy64Xbx4Uf7+/rp48aKrSwEAALDLHRfcAAAAzIrgBgAAYBIENwAAAJMguAEAAJgEwQ0AAMAkCG4AAAAmQXADAAAwCYIbAACASRDcAAAATILgBgAAYBIENwAAAJNwaXDbsGGDunfvrkqVKslisWjp0qU37L948WJ16tRJ5cqVk5+fn1q1aqUVK1YUUrUAAACu5dLgdvnyZTVu3FjvvfeeTf03bNigTp06KSoqStHR0brvvvvUvXt37d6928mVAgAAuJ7FMAzD1UVIksVi0ZIlS9SzZ0+7jmvQoIHCw8M1ceJEm/qnpKTI399fycnJ8vPzK0ipAAAALlHM1QXciqysLF28eFFlypTJt09aWprS0tKs2ykpKYVRGgC4XmqK5OEjuXu4uhIADmLqyQnTp0/X5cuX1bt373z7TJs2Tf7+/tZPUFBQIVYIAC5wZJ30fhvp1SDptRrS8vFSRtpNDwNQ9Jk2uM2fP1+TJ0/WwoULVb58+Xz7jR8/XsnJydbPiRMnCrFKAChkp2OlL3tLp/Zkb6dflLbNlqLGuLYuAA5hyuC2cOFCDRo0SF9//bU6dux4w75eXl7y8/PL8QGA29bPc6TMPEbXYuZLl88Vfj0AHMp0wW3+/Pnq37+/vvrqK3Xr1s3V5QBA0XLheN7tWdeklJOFWwsAh3NpcLt06ZJiYmIUExMjSTp69KhiYmIUHx8vKfs2Z9++fa3958+fr759+2r69Om65557dOrUKZ06dUrJyckuqR8AipxKTfJu9/KXyt5VuLUAcDiXBredO3eqadOmatq0qSRp1KhRatq0qXVpj4SEBGuIk6QPPvhAGRkZGjp0qCpWrGj9PP/88y6pHwDslpUpndghnYyWnLEa093PSL6BudvbjJQ8Szj+egAKVZFZx62wsI4bAJc5vFb679A/b1mWuUt6ZK5Uqaljr5MUL218Szq2SSpRTvrHIKnhI469BgCXILgBQGG4nCjNaCRdu5yzvUR5acReycPbNXUBMBXTTU4AAFP6dVHu0CZJl89IB3nnMgDbENwAoDBcvVCwfQDwFwQ3ACgMNe/Lu93iJtVsX5iVADAxghsAFIaqLaUmT+Zubz1CKl290MsBYE5MTgCAwmIY0m/LpLjvJbdiUsNe0l33u7oqACZCcAMAADAJbpUCAACYBMENAADAJAhuAAAAJkFwAwAAMAmCGwAAgEkQ3AAAAEyC4AYAAGASBDcAAACTILgBAACYBMENAADAJAhuAAAAJkFwAwAAMAmCGwAAgEkQ3AAAAEyC4AYAAGASBDcAAACTILgBAACYBMENAADAJAhuAAAAJkFwAwAAMAmCGwAAgEkQ3AAAAEyC4AYAAGASBDcAAACTILgBAACYBMENAADAJAhuAAAAJkFwAwAAMAmCGwAAgEkQ3AAAAEyC4AYAAGASBDcAAACTILgBAACYBMENAADAJAhuAAAAJkFwAwAAMAmCGwAAgEkQ3AAAAEyC4AYAAGASBDcAAACTILgBAACYBMENAADAJAhuAAAAJkFwAwAAMAmCGwAAgEkQ3AAAAEyC4AYAAGASLg1uGzZsUPfu3VWpUiVZLBYtXbr0psesX79ezZs3l7e3t2rWrKn333+/ECoFAABwPZcGt8uXL6tx48Z67733bOp/9OhRPfDAA2rTpo12796tf/3rXxo+fLgWLVrk5EoBAABcz2IYhuHqIiTJYrFoyZIl6tmzZ759XnjhBX333XeKi4uztkVEROiXX37R1q1bbbpOSkqK/P39lZycLD8/v1uuGwAAoLCY6hm3rVu3KiwsLEdb586dtXPnTl27ds1FVQEAABSOYq4uwB6nTp1SYGBgjrbAwEBlZGQoMTFRFStWzHVMWlqa0tLSrNspKSlOrxMAAMAZTDXiJmXfUv2r63d6/95+3bRp0+Tv72/9BAUFOb1GAAAAZzBVcKtQoYJOnTqVo+3MmTMqVqyYypYtm+cx48ePV3JysvVz4sSJwigVAADA4Ux1q7RVq1b6/vvvc7StXLlSLVq0kIeHR57HeHl5ycvLqzDKAwAAcCqXjrhdunRJMTExiomJkZS93EdMTIzi4+MlZY+W9e3b19o/IiJCx48f16hRoxQXF6d58+Zp7ty5GjNmjEvqBwAAKEwuHXHbuXOn7rvvPuv2qFGjJEn9+vXTJ598ooSEBGuIk6QaNWooKipKI0eO1KxZs1SpUiXNnDlTvXr1KvTaAQAACluRWcetsLCOGwAAMCtTTU4AAAC4kxHcAAAATILgBgAAYBIENwAAAJMguAEAAJgEwQ0AAMAkCG4AAAAmQXADAAAwCYIbAACASRDcAAAATILgBgAAYBIENwAAAJMguAEAAJgEwQ0AAMAkCG4AAAAmQXADAAAwCYIbAACASRDcAAAATILgBgAAYBIENwAAAJMguAEAAJgEwQ0AAMAkCG4AAAAmQXADAAAwCYIbAACASRDcAAAATILgBgAAYBIENwAAAJMguAEAAJgEwQ0AAMAkCG4AAAAmQXADAAAwCYIbAACASRDcAAAATILgBgAAYBIENwAAAJMoZk/n5ORkLVmyRBs3btSxY8d05coVlStXTk2bNlXnzp0VGhrqrDoBAADueDaNuCUkJOjpp59WxYoVNXXqVF2+fFlNmjRRhw4dVKVKFa1du1adOnVS/fr1tXDhQmfXDAAAcEeyacStcePG6tu3r3bs2KGQkJA8+1y9elVLly7VW2+9pRMnTmjMmDEOLRQAAOBOZzEMw7hZp7Nnz6pcuXI2n9Te/oUpJSVF/v7+Sk5Olp+fn6vLAQAAsJlNt0rtDWFFNbQBgGlcPC2d+lXKSHN1JQCKEJsnJ2zYsMGmfm3bti1wMQBwx0tNkb4bJsV9JxlZUvEAqcNEqXk/V1cGoAiw6VapJLm5uclisUiS8jvEYrEoMzPTcdU5AbdKARRp3w6Sfv02d3u/76Ua/MUYuNPZPOJWunRplSxZUv3791efPn0UEBDgzLoA4M5z5by0b0ne+3bOI7gBsH0B3oSEBL322mvaunWrGjZsqEGDBmnLli3y8/OTv7+/9QMAKKAr5yUjn7sWl84Wbi0AiiSbg5unp6fCw8O1YsUK7d+/X40aNdJzzz2noKAgTZgwQRkZGc6sEwBuf2VqSH6V895X/d7CrQVAkWTzM255OXr0qAYNGqT169fr7NmzKlOmjCNrcwqecQNQpP26SFr0VPbEhOvK3CU99ZNUvOj/PxaAc9n1yitJSktL06JFizRv3jxt3bpV3bp107Jly0wR2gCgyAvpJZWqnv1M28UEqVorqcUgQhsASXYEtx07dujjjz/WggULVKNGDfXv319ff/01gQ0AHK1K8+wPAPyNXcuBVK1aVf369VPz5vn/D+XBBx90WHHOwK1SAABgVnYFt5uejHXcAAAAnMbmW6VZWVk37wQAAACnsXk5kOvS0tJ0+fJlZ9QCAACAG7A5uCUmJqpbt27y9fWVn5+fQkNDdeTIEWfWBgAAgL+wObiNHz9e0dHRmjJlit544w0lJibq2WefveUCZs+erRo1asjb21vNmzfXxo0bb9j/yy+/VOPGjVW8eHFVrFhRAwYM0Llz5265DgAAgKLO5skJVatW1fvvv68HHnhAkvTbb78pJCREV69elYeHR4EuvnDhQvXp00ezZ89W69at9cEHH+ijjz5SbGysqlatmqv/pk2b1K5dO7399tvq3r27Tp48qYiICNWuXVtLluTzfr+/YXICAAAwK5uDW7FixXTixAlVrFjR2la8eHHFxcWpWrVqBbp4y5Yt1axZM0VGRlrbgoOD1bNnT02bNi1X/zfffFORkZE6fPiwte3dd9/V66+/rhMnTth0TYIbAAAwK5tvlRqGoWLFck5CLVasWIFnm6anpys6OlphYWE52sPCwrRly5Y8jwkNDdXvv/+uqKgoGYah06dP69tvv1W3bt0KVAMAAICZ2LwciGEY6tChQ47wduXKFXXv3l2enp7Wtl27dtl0vsTERGVmZiowMDBHe2BgoE6dOpXnMaGhofryyy8VHh6u1NRUZWRk6MEHH9S7776b73XS0tKUlpZm3U5JSbGpPgAAgKLG5uA2adKkXG09evS45QIsFkuObcMwcrVdFxsbq+HDh2vixInq3LmzEhISNHbsWEVERGju3Ll5HjNt2jRNmTLllusEAABwNZufcYuPj1eVKlVseoOCLdLT01W8eHF98803euihh6ztzz//vGJiYrR+/fpcx/Tp00epqan65ptvrG2bNm1SmzZt9Mcff+R4/u66vEbcgoKCeMYNAACYjs0prEaNGkpMTHTYhT09PdW8eXOtWrUqR/uqVasUGhqa5zFXrlzJFRzd3d0lZY/U5cXLy0t+fn45PgAAAGZk1+QERxs1apQ++ugjzZs3T3FxcRo5cqTi4+MVEREhKXvtuL59+1r7d+/eXYsXL1ZkZKSOHDmizZs3a/jw4br77rtVqVIlh9cHAABQlNj8jJszhIeH69y5c5o6daoSEhIUEhKiqKgo6/IiCQkJio+Pt/bv37+/Ll68qPfee0+jR49WqVKldP/99+u1115z1VcAAAAoNDY/4+bm5qZXXnlFvr6+N+w3fPhwhxTmLKzjBgAAzMqu4FalShXrM2V5nsxiKfLvLyW4AQAAs7LrVunOnTtVvnx5Z9UCAACAG7B5ckJ+a6sBAACgcLh0VikAAABsZ3NwmzRp0k0nJgAAAMB5bJqcEB8fr6pVq9p80pMnT6py5cq3VJizMDkBAACYlU0jbv/4xz/0zDPPaMeOHfn2SU5O1pw5cxQSEqLFixc7rEAAwC24lirt+kxa/Ky08t/S2QOurgjALbBpVmlcXJz+85//qEuXLvLw8FCLFi1UqVIleXt768KFC4qNjdW+ffvUokULvfHGG+ratauz6wYA3Ez6ZenT7tLJ6D/btr0v9f5MqveA6+oCUGA2r+MmSampqYqKitLGjRt17NgxXb16VQEBAWratKk6d+6skJAQZ9bqENwqBXDH2DpbWjE+d7t/kPT8HsnN5secARQRdgW32wHBDcAd4/OHpcOr8943eKsUWL9w6wFwy/jrFgDcrrz9C7YPQJFFcAOA21WzPnm333W/5F80Z/4DuDGCGwDcru66Xwp7RfIo8Wdb1VZSz/ddVxOAW8IzbgBwu0tNlk7uknzLS4ENXF0NgFtg10vmAQAm5O0v3XXfjftkpEvFPAunHgAFVqDgduDAAa1bt05nzpxRVlZWjn0TJ050SGEAgELw81xp0wwpOV4qW1tq/6LU8BFXVwUgH3bfKp0zZ44GDx6sgIAAVahQQRaL5c+TWSzatWuXw4t0JG6VAsD/RH8iff987vbwL6XgfxZ6OQBuzu7gVq1aNQ0ZMkQvvPCCs2pyKoIbAPzPzKbS+SO524PukQatKPx6ANyU3bNKL1y4oEcffdQZtQAACothSOeP5r3v/OHCrQWAzewObo8++qhWrlzpjFoAAIXFYpEq5POawgqNCrcWADaze3JCrVq19NJLL2nbtm1q2LChPDw8cuwfPny4w4oDADhR+/HSgicl/eWJGTcPqe0Yl5UE4MbsfsatRo0a+Z/MYtGRI3k8L1GE8IwbAPzFoZ+kzTOlc4ez13hrM0qqeo+rqwKQDxbgBQAAMIlbeuWVYRi6w3IfAACAyxQouH322Wdq2LChfHx85OPjo0aNGunzzz93dG0AAAD4C7snJ7z11lt66aWX9Nxzz6l169YyDEObN29WRESEEhMTNXLkSGfUCQAAcMcr0OSEKVOmqG/fvjnaP/30U02ePFlHj+azLlARwTNuAG5bSfHSyWjJP0iq0sLV1QBwArtH3BISEhQaGpqrPTQ0VAkJCQ4pCgBgB8OQfhwn/fyRZPzv/dGVW0iPL5B8y7m2NgAOZfczbrVq1dLXX3+dq33hwoWqXbu2Q4oCANgh5ktpx4d/hjZJOrlT+mGE62oC4BR2j7hNmTJF4eHh2rBhg1q3bi2LxaJNmzZp9erVeQY6AICT/bIg7/b9P0qpyZK3f+HWA8Bp7B5x69Wrl7Zv366AgAAtXbpUixcvVkBAgHbs2KGHHnrIGTUCAG7k2pW8241MKSOtcGsB4FR2j7hJUvPmzfXFF184uhYAQEHU6ZI9KeHvKjWTfMsXfj0AnMam4JaSkmKdgZmSknLDvszUBIBCds9gaX+U9MfuP9u8/aVub7quJgBOYdNyIO7u7kpISFD58uXl5uYmi8WSq49hGLJYLMrMzHRKoY7CciAAbksZ6VLsf6Xff5b8K0uNn2BGKXAbsmnEbc2aNSpTpowkae3atU4tCABQAMU8pUaPZn8A3LZ4yTwAAIBJ2D2rdPny5dq0aZN1e9asWWrSpImeeOIJXbhwwaHFAQAA4E92B7exY8daJyjs3btXo0aN0gMPPKAjR45o1KhRDi8QAAAA2exeDuTo0aOqX7++JGnRokXq3r27/vOf/2jXrl164IEHHF4gAAAAstk94ubp6akrV7IXe/zpp58UFhYmSSpTpsxNlwoBAABAwdk94nbvvfdq1KhRat26tXbs2KGFCxdKkg4cOKAqVao4vEAAAABks3vE7b333lOxYsX07bffKjIyUpUrV5Yk/fjjj+rSpYvDCwQAAEA2lgMBAAAwCbtH3Hbt2qW9e/dat//73/+qZ8+e+te//qX09HSHFgcAAIA/2R3cnn32WR04cECSdOTIET322GMqXry4vvnmG40bN87hBQIAACCb3cHtwIEDatKkiSTpm2++Udu2bfXVV1/pk08+0aJFixxeIAAAALLZHdwMw1BWVpak7OVArq/dFhQUpMTERMdWBwAAACu7g1uLFi30yiuv6PPPP9f69evVrVs3SdkL8wYGBjq8QAAAAGSzO7jNmDFDu3bt0nPPPacJEyaoVq1akqRvv/1WoaGhDi8QAAAA2Ry2HEhqaqrc3d3l4eHhiNM5DcuBAAAAs7J7xE2SkpKS9NFHH2n8+PE6f/68JCk2NlZnzpxxaHEAgJu4dFbaMUfa8p507rCrqwHgZHaPuO3Zs0cdOnRQqVKldOzYMe3fv181a9bUSy+9pOPHj+uzzz5zVq0OwYgbgNtG3A/SokFSRur/GizSfROkdmNdWhYA57F7xG3UqFEaMGCADh48KG9vb2t7165dtWHDBocWBwDIR9olaUnEX0KbJBnS2lekhD0uKwuAc9kd3H7++Wc9++yzudorV66sU6dOOaQoAMBNHF4tpV/Me1/s0sKtBUChsTu4eXt7KyUlJVf7/v37Va5cObsLmD17tmrUqCFvb281b95cGzduvGH/tLQ0TZgwQdWqVZOXl5fuuusuzZs3z+7rAoC5WQq4D3COdevWyWKxKCkpydWl3NbsDm49evTQ1KlTde3aNUmSxWJRfHy8XnzxRfXq1cuucy1cuFAjRozQhAkTtHv3brVp00Zdu3ZVfHx8vsf07t1bq1ev1ty5c7V//37Nnz9f9erVs/drAIC51eogefnnva/BQ4VbCyApNDRUCQkJ8vfP579LOITdkxNSUlL0wAMPaN++fbp48aIqVaqkU6dOqVWrVoqKilKJEiVsPlfLli3VrFkzRUZGWtuCg4PVs2dPTZs2LVf/5cuX67HHHtORI0dUpkwZe8rOUT+TEwDcFvb/KH07ULp2JXvb4iZ1mCjdO9K1dQFwGrtH3Pz8/LRp0yYtWrRIr776qp577jlFRUVp/fr1doW29PR0RUdHKywsLEd7WFiYtmzZkucx3333nVq0aKHXX39dlStXVp06dTRmzBhdvXrV3q8BAOZXt6s0cp/UfabU9XVp+G5CGxymffv2GjZsmEaMGKHSpUsrMDBQH374oS5fvqwBAwaoZMmSuuuuu/Tjjz9Kyn2r9Pjx4+revbtKly6tEiVKqEGDBoqKirKePzY2Vg888IB8fX0VGBioPn36WF+duW7dOnl6euZ4fGr69OkKCAhQQkKCtb7nnntOzz33nEqVKqWyZcvq3//+txy0PG2RVcyezhkZGfL29lZMTIzuv/9+3X///QW+cGJiojIzM3O9JiswMDDfSQ5HjhzRpk2b5O3trSVLligxMVFDhgzR+fPn833OLS0tTWlpadbtvJ7PAwDTKl5Gat7P1VXgNvXpp59q3Lhx2rFjhxYuXKjBgwdr6dKleuihh/Svf/1Lb7/9tvr06ZPnI05Dhw5Venq6NmzYoBIlSig2Nla+vr6SpISEBLVr105PP/203nrrLV29elUvvPCCevfurTVr1qh9+/YaMWKE+vTpo19++UXHjh3ThAkTNH/+fFWsWDFHfYMGDdL27du1c+dOPfPMM6pWrZqefvrpQvsdFTrDTjVr1jRiYmLsPSyXkydPGpKMLVu25Gh/5ZVXjLp16+Z5TKdOnQxvb28jKSnJ2rZo0SLDYrEYV65cyfOYSZMmGZJyfZKTk2/5OwAAcLtq166dce+991q3MzIyjBIlShh9+vSxtiUkJBiSjK1btxpr1641JBkXLlwwDMMwGjZsaEyePDnPc7/00ktGWFhYjrYTJ04Ykoz9+/cbhmEYaWlpRtOmTY3evXsbDRo0MJ566qlc9QUHBxtZWVnWthdeeMEIDg6+tS9exNl9q/Tf//53jjcmFFRAQIDc3d1zja6dOXMm35fVV6xYUZUrV87x4GNwcLAMw9Dvv/+e5zHjx49XcnKy9XPixIlbqhsAgDtFo0aNrD+7u7urbNmyatiwobXt+p/Xeb05afjw4XrllVfUunVrTZo0SXv2/Lm+YHR0tNauXStfX1/r5/pEw8OHs98A4unpqS+++EKLFi3S1atXNWPGjFzXuOeee2Sx/DmLulWrVjp48KAyMzNv8ZsXXXYHt5kzZ2rjxo2qVKmS6tatq2bNmuX42MrT01PNmzfXqlWrcrSvWrUq35fVt27dWn/88YcuXbpkbTtw4IDc3NxUpUqVPI/x8vKSn59fjg8AALi5v79/3GKx5Gi7HpqysrJyHfvUU0/pyJEj6tOnj/bu3asWLVro3Xfftfbv3r27YmJicnwOHjyotm3bWs9x/Zn38+fP3/KA0e3CrmfcJKlnz54Ou/ioUaPUp08ftWjRQq1atdKHH36o+Ph4RURESMoeLTt58qT1NVpPPPGEXn75ZQ0YMEBTpkxRYmKixo4dq4EDB8rHx8dhdQEAgFsXFBSkiIgIRUREaPz48ZozZ46GDRumZs2aadGiRapevbqKFcs7ihw+fFgjR47UnDlz9PXXX6tv375avXq13Nz+HHPatm1bjmO2bdum2rVry93d3anfy5XsDm6TJk1y2MXDw8N17tw5TZ06VQkJCQoJCVFUVJSqVasmKfvhxb8+8Ojr66tVq1Zp2LBhatGihcqWLavevXvrlVdecVhNAADg1o0YMUJdu3ZVnTp1dOHCBa1Zs0bBwcGSsicuzJkzR48//rjGjh2rgIAAHTp0SAsWLNCcOXMkSX369FFYWJgGDBigrl27qmHDhpo+fbrGjv3zXbwnTpzQqFGj9Oyzz2rXrl169913NX36dJd838Jid3C7bufOnYqLi5PFYlFwcLCaN29eoPMMGTJEQ4YMyXPfJ598kqutXr16uW6vAgCAoiUzM1NDhw7V77//Lj8/P3Xp0kVvv/22JKlSpUravHmzXnjhBXXu3FlpaWmqVq2aunTpIjc3N7388ss6duyYvv/+e0lShQoV9NFHH6l3797q1KmTmjRpIknq27evrl69qrvvvlvu7u4aNmyYnnnmGZd958Jg9wK8v//+ux5//HFt3rxZpUqVkiQlJSUpNDRU8+fPV1BQkFMKdRQW4AUAwPzat2+vJk2a5Dlp4XZm9+SEgQMH6tq1a4qLi7M+LBgXFyfDMDRo0CBn1AgAAAAV4Fbpxo0btWXLFtWtW9faVrduXb377rtq3bq1Q4sDAADAn+wOblWrVrW+YP6vMjIyVLlyZYcUBQAAcCPr1q1zdQkuYfet0tdff13Dhg3Tzp07re8D27lzp55//nm9+eabDi8QAAAA2eyenFC6dGlduXJFGRkZ1rVXrv/895fMF8XF8picAAAAzMruW6V32uwNAACAosLuETezY8QNAACYld3PuM2dOzfP9oyMDI0fP/6WCwIAAEDe7A5uo0ePVq9evXI8v/bbb7/p7rvv1tdff+3Q4gAAAPAnu4Pb7t27dfr0aTVs2FCrVq3SrFmz1KxZM4WEhCgmJsYZNQIAAORr3bp1slgsSkpKcnUpTmd3cKtRo4Y2bNigRx55RF26dNHIkSM1b948ffbZZypZsqQzagQAAIXo1KlTGjZsmGrWrCkvLy8FBQWpe/fuWr16tcOu0b59e40YMcJh57tTFOgl8z/88IPmz5+v0NBQ7d+/X3PmzFHbtm1VqVIlR9cHAMAdKzPL0I6j53XmYqrKl/TW3TXKyN3N4tRrHjt2TK1bt1apUqX0+uuvq1GjRrp27ZpWrFihoUOH6rfffnPq9f/KMAxlZmZalx9DAUbcnn32WfXu3Vvjxo3Thg0btGfPHnl5ealhw4Y84wYAgIMs/zVB9762Ro/P2abnF8To8TnbdO9ra7T81wSnXnfIkCGyWCzasWOHHnnkEdWpU0cNGjTQqFGjtG3bNklSfHy8evToIV9fX/n5+al37946ffq09RyTJ09WkyZN9Pnnn6t69ery9/fXY489posXL0qS+vfvr/Xr1+udd96RxWKRxWLRsWPHrLc8V6xYoXXt5IgAACAASURBVBYtWsjLy0sbN25UWlqahg8frvLly8vb21v33nuvfv75Z6f+Hooqu4Pb5s2btX37do0ZM0YWi0UVKlRQVFSUpk6dqoEDBzqjRgAA7ijLf03Q4C92KSE5NUf7qeRUDf5il9PC2/nz57V8+XINHTo016L6klSqVCkZhqGePXvq/PnzWr9+vVatWqXDhw8rPDw8R9/Dhw9r6dKl+uGHH/TDDz9o/fr1evXVVyVJ77zzjlq1aqWnn35aCQkJSkhIUFBQkPXYcePGadq0aYqLi1OjRo00btw4LVq0SJ9++ql27dqlWrVqqXPnzkVyoX9nszu4RUdHq3Hjxrnahw4dqujoaIcUBQDAnSozy9CU72OV1yKr19umfB+rzCzHL8N66NAhGYahevXq5dvnp59+0p49e/TVV1+pefPmatmypT7//HOtX78+xyhYVlaWPvnkE4WEhKhNmzbq06eP9Rk5f39/eXp6qnjx4qpQoYIqVKggd3d367FTp05Vp06ddNddd8nb21uRkZF644031LVrV9WvX19z5syRj49PvkuU3c5sDm5nzpyRJHl5eeW5PyMjQ8nJyY6pCgCAO9SOo+dzjbT9lSEpITlVO446frTp+pr8Fkv+z9HFxcUpKCgoxwhZ/fr1VapUKcXFxVnbqlevnmPSYsWKFa1Z4mZatGhh/fnw4cO6du2aWrdubW3z8PDQ3XffneN6dwqbg9vff+HBwcGKj4+3bp87d06tWrVybHUAANxhzlzMP7QVpJ89ateuLYvFcsNAZBhGnsHu7+0eHh459lssFmVlZdlUx19v0+YXJvOr43Znc3D7+5uxfv/9d2VkZNywDwAAsE/5kt4O7WePMmXKqHPnzpo1a5YuX76ca39SUpLq16+v+Ph4nThxwtoeGxur5ORkBQcH23wtT09PZWZm3rRfrVq15OnpqU2bNlnbrl27pp07d9p1vduF3c+43cidmHwBAHCku2uUUUV/b+X3J6pFUkX/7KVBnGH27NnKzMzU3XffrUWLFungwYOKi4vTzJkz1apVK3Xs2FGNGjXSk08+qV27dmnHjh3q27ev2rVrl+MW581Ur15d27dv17Fjx5SYmJjvaFyJEiU0ePBgjR07VsuXL1dsbKyefvppXblyRYMGDXLU1zYNhwY3AABwa9zdLJrUvb4k5Qpv17cnda/vtPXcatSooV27dum+++7T6NGjFRISok6dOmn16tWKjIyUxWLR0qVLVbp0abVt21YdO3ZUzZo1tXDhQruuM2bMGLm7u6t+/foqV65cjsev/u7VV19Vr1691KdPHzVr1kyHDh3SihUrVLp06Vv9uqZjMWy8v+nu7q4DBw6oXLlyMgxDQUFB2rRpk6pXry5JOn36tOrVq2fTsKcrpaSkyN/fX8nJyfLz83N1OQAA5Gn5rwma8n1sjokKFf29Nal7fXUJqejCyuBKNgc3Nze3HLdC//5Q4PVtghsAAI7hijcnoGiz+R0Sa9eudWYdAADgb9zdLGp1V1lXl4EixOYRt9sFI24AAMCsmJwAAABgEgQ3AAAAkyC4AQAAmATBDQAAwCQIbgAAACZhV3Bbu3atpk+frs2bN0uSPvjgA1WtWlXlypXT008/ratXrzqlSAAAANgR3ObMmaNOnTopMjJSHTp00LRp0zR69Gh169ZNvXv31tdff60pU6Y4s1YAAFCEtG/fXiNGjHB1GabTv39/9ezZs0DH2hzc3nnnHb399ts6dOiQli5dqokTJ2rWrFmKjIzUrFmz9NFHH+nbb78tUBEAAKBo6N+/vywWiyIiInLtGzJkiCwWi/r37y9JWrx4sV5++eVCrvDOZnNwO3LkiB588EFJUpcuXWSxWHT33Xdb97ds2VInTpxwfIUAANypsjKloxulvd9m/zOrcF4rGRQUpAULFuR4BCo1NVXz589X1apVrW1lypRRyZIlC6WmgsjMzFRWVpary3Aom4NbamqqfHx8rNteXl7y8vLKsZ2RkeHY6gAAuFPFfifNCJE+/ae0aFD2P2eEZLc7WbNmzVS1alUtXrzY2rZ48WIFBQWpadOm1ra/3yqdPXu2ateuLW9vbwUGBuqRRx6x7jMMQ6+//rpq1qwpHx8fNW7c2HqnzjAMdezYUV26dNH1FzolJSWpatWqmjBhgiRp3bp1slgsWrZsmRo3bixvb2+1bNlSe/futV7jk08+UalSpfTDDz+ofv368vLy0vHjx5Wenq5x48apcuXKKlGihFq2bKl169ZZjzt+/Li6d++u0qVLq0SJEmrQoIGioqIkSRcuXNCTTz6pcuXKycfHR7Vr19bHH39sPfbkyZMKDw9X6dKlVbZsWfXo0UPHjh2z7s/MzNSoUaNUqlQplS1bVuPGjdOtvLTK5uBmsVh08eJFpaSkKDk5WRaLRZcuXVJKSor1AwAAHCD2O+nrvlLKHznbUxKy2wshvA0YMCBHQJk3b54GDhyYb/+dO3dq+PDhmjp1qvbv36/ly5erbdu21v3//ve/9fHHHysyMlL79u3TyJEj9X//939av369LBaLPv30U+3YsUMzZ86UJEVERCgwMFCTJ0/OcZ2xY8fqzTff1M8//6zy5cvrwQcf1LVr16z7r1y5omnTpumjjz7Svn37VL58eQ0YMECbN2/WggULtGfPHj366KPq0qWLDh48KEkaOnSo0tLStGHDBu3du1evvfaafH19JUkvvfSSYmNj9eOPPyouLk6RkZEKCAiwXuu+++6Tr6+vNmzYoE2bNsnX11ddunRRenq6JGn69OmaN2+e5s6dq02bNun8+fNasmRJwf/FGDayWCyGm5ub9ZPfdlGXnJxsSDKSk5NdXQoAALllZhjG9HqGMckvn4+/YUwPzu7nBP369TN69OhhnD171vDy8jKOHj1qHDt2zPD29jbOnj1r9OjRw+jXr59hGIbRrl074/nnnzcMwzAWLVpk+Pn5GSkpKbnOeenSJcPb29vYsmVLjvZBgwYZjz/+uHX766+/Nry8vIzx48cbxYsXN/bv32/dt3btWkOSsWDBAmvbuXPnDB8fH2PhwoWGYRjGxx9/bEgyYmJirH0OHTpkWCwW4+TJkzmu3aFDB2P8+PGGYRhGw4YNjcmTJ+f5++jevbsxYMCAPPfNnTvXqFu3rpGVlWVtS0tLM3x8fIwVK1YYhmEYFStWNF599VXr/mvXrhlVqlQxevTokec5b6aYrQFv7dq1BU+HAADANse35B5py8GQUk5m96vRxmllBAQEqFu3bvr0009lGIa6detmHWnKS6dOnVStWjXVrFlTXbp0UZcuXfTQQw+pePHiio2NVWpqqjp16pTjmPT09By3Xh999FEtWbJE06ZNU2RkpOrUqZPrOq1atbL+XKZMGdWtW1dxcXHWNk9PTzVq1Mi6vWvXLhmGketcaWlpKlu2rCRp+PDhGjx4sFauXKmOHTuqV69e1nMMHjxYvXr10q5duxQWFqaePXsqNDRUkhQdHa1Dhw7les4vNTVVhw8fVnJyshISEnLUXKxYMbVo0aLAt0ttDm7t2rUr0AUAAIAdLp12bL9bMHDgQD333HOSpFmzZt2wb8mSJbVr1y6tW7dOK1eu1MSJEzV58mT9/PPP1gkCy5YtU+XKlXMc99fn5a9cuaLo6Gi5u7tbb2PawmKxWH/28fHJsZ2VlSV3d3fref/q+u3Qp556Sp07d9ayZcu0cuVKTZs2TdOnT9ewYcPUtWtXHT9+XMuWLdNPP/2kDh06aOjQoXrzzTeVlZWl5s2b68svv8xVU7ly5Wyu3x4Oe3NCRkaG4uPjHXU6AADuTL6Bju13C64/q5Wenq7OnTvftH+xYsXUsWNHvf7669qzZ4+OHTumNWvWWCcKxMfHq1atWjk+QUFB1uNHjx4tNzc3/fjjj5o5c6bWrFmT6xrbtm2z/nzhwgUdOHBA9erVy7empk2bKjMzU2fOnMl17QoVKlj7BQUFKSIiQosXL9bo0aM1Z84c675y5cqpf//++uKLLzRjxgx9+OGHkrIncRw8eFDly5fPdW5/f3/5+/urYsWKOWrOyMhQdHT0TX+X+bF5xO1m9u3bp2bNmikzs3CmKgMAcFuqFir5VcqeiKC8bqdZsvdXC3V6Ke7u7tbbkH8frfq7H374QUeOHFHbtm1VunRpRUVFKSsrS3Xr1lXJkiU1ZswYjRw5UllZWbr33nuVkpKiLVu2yNfXV/369dOyZcs0b948bd26Vc2aNdOLL76ofv36ac+ePSpdurT1OlOnTlXZsmUVGBioCRMmKCAg4IaL2dapU0dPPvmk+vbtq+nTp6tp06ZKTEzUmjVr1LBhQz3wwAMaMWKEunbtqjp16ujChQtas2aNgoODJUkTJ05U8+bN1aBBA6WlpemHH36w7nvyySf1xhtvqEePHpo6daqqVKmi+Ph4LV68WGPHjlWVKlX0/PPP69VXX1Xt2rUVHByst956S0lJSQX+d8K7SgEAyswq+PIEcDA3d6nLa//bsPxt5/+2u7ya3a8Q+Pn5yc/P76b9SpUqpcWLF+v+++9XcHCw3n//fc2fP18NGjSQJL388suaOHGipk2bpuDgYHXu3Fnff/+9atSoobNnz2rQoEGaPHmymjVrJkmaNGmSKlWqlGsh4FdffVXPP/+8mjdvroSEBH333Xfy9PS8YW0ff/yx+vbtq9GjR6tu3bp68MEHtX37dutoX2ZmpoYOHarg4GB16dJFdevW1ezZsyVlPzM3fvx4NWrUSG3btpW7u7sWLFggSSpevLg2bNigqlWr6uGHH1ZwcLAGDhyoq1evWn9no0ePVt++fdW/f3+1atVKJUuW1EMPPWTHv4GcLIaNT8dd/0Xm5+rVqzpw4ECRH3FLSUmRv7+/kpOTbfoPEQBuZ9uOnNPry3/TrvgklS3hqf+7p5qGd6gtd7e/BwYUutjvpOUv5Jyo4Fc5O7TVf9B1dbnIunXrdN999+nChQsqVaqUq8txGZtvlcbGxuqxxx5TjRo18tyfkJCgAwcOOKwwAIBz/XYqRX3n7VB6RvaD4+cup+ud1Qd1KS1DL/2zvourg+o/KNXrlj179NLp7GfaqoUW2kgbiiabg1tISIhatmypwYMH57k/JiYmx4N8AICi7ZPNx6yh7a++3H5cIzrWVklvDxdUhRzc3J265AfMx+bgdu+992r//v357i9ZsmSOFZIBALaJ/SNFs9YeUsyJJFUp7aNB99ZQWIMKNz/wFh1NvJxne+q1LJ1OSSO4oUhp3779Lb0q6nZh8zNutwuecQNQlPx2KkUPz96iK+k5nw9+89HGeqR5Fadee8r3+/Tx5mO52v28i2nHhI7y9uCWHFDUMKsUAFwoct3hXKFNkmb8dMDpowsDW9eQv0/uUbWI9ncR2oAiyuZbpbYurlu1atUCFwMAd5p9f6Tk2f77hatKunJNpUvceJmDWxFUprgWDW6lmasPacfR8yrv56U+91TToy2Cbn4wAJewObhVr149xyskrjMMw9pusViUkZHhuOoA4DZXrUxxHTpzKVd7mRKeKuntsDXS81WrfEnNfLzpzTsCKBJs/r/C7t2782w3DEMLFizQzJkzre/8AgDYZlCbGlq7/4z+vv7twNbVVcydp1kA5HRLkxN++uknvfjiizpw4IBGjRqlMWPGFPnwxuQEAEXNyn2nNH3lAe0/fVHlSnppYOsaimhXM8+7HADubAUah4+OjtaLL76ojRs36qmnnlJUVJTKly/v6NoA4I4Q1qCCwhpUUOq1THkVcyOwAciXXePwhw4dUnh4uFq2bKly5copNjZW7733HqENABzA28Od0AbghmwObkOGDFGDBg2UnJysnTt36quvvlLNmjWdWRsAAAD+wuZn3Nzc3OTt7a169erdsN+uXbvsKmD27Nl64403lJCQoAYNGmjGjBlq0+bmr/fYvHmz2rVrp5CQEMXExNh8PZ5xAwAAZmXzM26TJk1y+MUXLlyoESNGaPbs2WrdurU++OADde3aVbGxsTdcDy45OVl9+/ZVhw4ddPr0aYfXBQAAUBS59JVXLVu2VLNmzRQZGWltCw4OVs+ePTVt2rR8j3vsscdUu3Ztubu7a+nSpYy4AQCAO4JdkxO2b9+uCRMmaNy4cVq5cuUtXTg9PV3R0dEKCwvL0R4WFqYtW7bke9zHH3+sw4cP2zwCmJaWppSUlBwfAAAAM7I5uC1ZskStW7fWO++8ow8//FBdu3bVjBkzCnzhxMREZWZmKjAwMEd7YGCgTp06lecxBw8e1Isvvqgvv/xSxYrZdpd32rRp8vf3t36CgniVCwAAMCebg9t//vMf9e/fX0lJSUpKStKUKVP0yiuv3HIBf5/6/tdXaP1VZmamnnjiCU2ZMkV16tSx+fzjx49XcnKy9XPixIlbrhkAAMAVbH7Gzc/PTzt37rSGprS0NJUoUUKnTp1SQECA3RdOT09X8eLF9c033+ihhx6ytj///POKiYnR+vXrc/RPSkpS6dKl5e7ubm3LysqSYRhyd3fXypUrdf/999/0ujzjBgAAzMrmEbdLly6pVKlS1m0vLy/5+PgU+JkxT09PNW/eXKtWrcrRvmrVKoWGhubq7+fnp7179yomJsb6iYiIUN26dRUTE6OWLVsWqA4AAACzsOuVVytWrJC/v791OysrS6tXr9avv/5qbXvwwQdtPt+oUaPUp08ftWjRQq1atdKHH36o+Ph4RURESMq+zXny5El99tlncnNzU0hISI7jy5cvL29v71ztAAAAtyO7glu/fv1ytT377LPWny0WizIzM20+X3h4uM6dO6epU6cqISFBISEhioqKUrVq1SRJCQkJio+Pt6dEAACA25ZL13FzBZ5xAwAAZmXXOm4AAABwHbuD2zfffKOHH35YISEhatiwoR5++GF9++23zqgNAAAAf2FzcMvKylJ4eLjCw8MVGxurWrVqqWbNmtq3b5/Cw8P12GOP6Q676woAAFCobJ6cMGPGDP3000/67rvv9M9//jPHvu+++04DBgzQO++8oxEjRji8SAAAANgxOaFRo0YaMWKEBg4cmOf+uXPnasaMGdq7d69DC3Q0JicAAACzsvlW6cGDB9WxY8d893fs2FGHDh1ySFEAAADIzebg5uPjo6SkpHz3p6SkyMfHxyFFAQAAIDebg1urVq0UGRmZ7/5Zs2apVatWDikKAAAAudk8OWHChAlq3769zp07pzFjxqhevXoyDENxcXGaPn26/vvf/2rt2rXOrBUAAOCOZtebE5YsWaJnnnlG58+fz9FeunRpffDBB+rVq5fDC3Q0JicAAACzsvuVV1euXNGKFSt08OBBSVKdOnUUFham4sWLO6VARyO4AQAAs3Lou0pPnjypypUrO+p0TkFwAwAAZuWQd5WeOnVKw4YNU61atRxxOgAAAOTB5uCWlJSkJ598UuXKlVOlSpU0c+ZMZWVlaeLEiapZs6a2bdumefPmObNWAACAO5rNs0r/9a9/acOGDerXr5+WL1+ukSNHavny5UpNTdWPP/6odu3aObNOAACAO57NwW3ZsmX6+OOP1bFjRw0ZMkS1atVSnTp1NGPGDGfWBwAAgP+x+VbpH3/8ofr160uSatasKW9vbz311FNOKwwAAAA52RzcsrKy5OHhYd12d3dXiRIlnFIUAAAAcrP5VqlhGOrfv7+8vLwkSampqYqIiMgV3hYvXuzYCgEAACDJjuDWr1+/HNv/93//5/BiAOB2djU9U+5uFnkWc8hKTADuQA5dgNcMWIAXQGE7ePqipnwfq82HE+Xh5qZujSpq4j/rq3QJT1eXBsBkbB5xAwDYL/nqNT0+Z7sSL6VJktIzs7Rk90nFn7+iRYNDXVwdALNhvB6mkZllKDPrjhogxm1g6e6T1tD2V9HHL2hX/AUXVATAzBhxQ5F39mKaXlkWqx/3nlKWYahzgwr69z+DVdHfx9WlATcVf/5KvvtOnL+iZlVLF2I1AMyO4IYiLSvLUJ+52/XbqYvWtmV7ExSbkKKVI9vKw51BYxRtDSv757uvQaX89wFAXvhTD0Xa+gNnc4S2644mXtbKfaddUBFgnwcaVlT9irknQj3UtLJqlfd1QUUAzIwRNxRpx85dLtA+oKjwLOam+U/fo8j1h/VT3Gl5e7ipZ5PK6h9a3dWlATAhghuKtLxGKqz7KrGcC8zBv7iHXuxaTy92refqUgCYHMENRc6nW47p063HdDYlTS2ql1bjKv765ffkHH2aVS2ldrXLuaZAAABchOCGIuWtVQc0c/VB6/ba/WdV3MNNfe6pqs2HzinLMNS1YUUNva+W3NwsLqwUAIDCR3BDkXE1PVPzNh3N1X7lWpYki9aMaV/oNQEAUJQwqxRFxsmkq7qUlpHnvv2nc88sBQDgTkNwQ5FRqZS3Sni657mvNssmAABAcEPRUdyzmAa0rpFHu3ue7QAA3Gl4xg1FyuiwOipV3EOfbzuu0ymp+kf1MhoTVpeFSgEAkGQxDOOOemt3SkqK/P39lZycLD8/1gEDAADmwa1SAAAAkyC4AQAAmATBDQAAwCQIbgAAACZBcAMAADAJghsAAIBJENwAAABMguCGW3aHLQUIAIDL8OYEFNjRxMv6f8vitG7/GXl7uOvhZpU1rks9+XrxnxUAAM7An7AokOSr1xT+wVaduZgmSbqUlqHPth7XsXNX9NnAu11cHQAAtydulaJAluz63Rra/mrDgbOK/SNF6RlZ2njwrDYdTNS1zCwXVAgAwO2HETcUyNHEy/nui9qboAU/xyvxUrokKcDXSzMfa6LQWgGFVR4AALclRtxQIMEV/fLd9+HGI9bQJkmJl9L07OfRSkm9VhilAQBw2yK4oUBCa5XNd196Ru5boxfTMrTi11POLAkAgNsewQ0Fcik10/5j0jKcUAkAAHcOghsK5K7yJVSmhKfN/S0WqX3d8k6sCACA2x/BDQXiVcxdY8Lq5mqvWqa4+reqnqt9aPtaqhFQohAqAwDg9sWsUhTYEy2rqmqZ4vpy+3GdvZime2qW1YDW1VXW10tdG1ZQ1N4EWSwWdWtUUf+oXsbV5QIAYHoWw8XvK5o9e7beeOMNJSQkqEGDBpoxY4batGmTZ9/FixcrMjJSMTExSktLU4MGDTR58mR17tzZ5uulpKTI399fycnJ8vPLf2YkAABAUePSW6ULFy7UiBEjNGHCBO3evVtt2rRR165dFR8fn2f/DRs2qFOnToqKilJ0dLTuu+8+de/eXbt37y7kygEAAAqfS0fcWrZsqWbNmikyMtLaFhwcrJ49e2ratGk2naNBgwYKDw/XxIkTberPiBsAADArl424paenKzo6WmFhYTnaw8LCtGXLFpvOkZWVpYsXL6pMmfyfn0pLS1NKSkqODwAAgBm5LLglJiYqMzNTgYGBOdoDAwN16pRtC7VOnz5dly9fVu/evfPtM23aNPn7+1s/QUFBt1Q3AACAq7h8ORCLxZJj2zCMXG15mT9/viZPnqyFCxeqfPn81wcbP368kpOTrZ8TJ07ccs0AAACu4LLlQAICAuTu7p5rdO3MmTO5RuH+buHChRo0aJC++eYbdezY8YZ9vby85OXldcv1AgAAuJrLRtw8PT3VvHlzrVq1Kkf7qlWrFBoamu9x8+fPV//+/fXVV1+pW7duzi4TAACgyHDpAryjRo1Snz591KJFC7Vq1Uoffvih4uPjFRERISn7NufJkyf12WefScoObX379tU777yje+65xzpa5+PjI39/f5d9DwAAgMLg0uAWHh6uc+fOaerUqUpISFBISIiioqJUrVo1SVJCQkKONd0++OADZWRkaOjQoRo6dKi1vV+/fvrkk08Ku3wAAIBC5fI3JxQ21nEDAABm5fJZpQAAALANwQ0AAMAkCG4AAAAmQXADAAAwCYIbAACASRDcAAAATMKl67jh9vf7hSvadDBRfj4eur9eeXl7uLu6JAAATIvgBqd5e9UBvbvmoLL+t1JggK+nPur3DzUJKuXawgAAMClulcIpth05p3dW/xnaJCnxUrqGzd+lrKw7as1nAAAchuAGp/julz/ybD9x/qp2n7hQyNUAAHB7ILjBKTIys/Lddy2TETcAAAqC4Aan6BJSIc/2AF8vNa9W2rp9LTNLcQkpOp2SWlilAQBgWkxOgFPcV7e8wlsEaeHOE9Y2bw83vfFoI3m4Z/99Yenuk/p/UXE6ezFNFovUMThQbz7SWP7FPQp0zYzMLCVdvaZSPh4q5s7fSQAAtx+LYRh31H2rlJQU+fv7Kzk5WX5+fq4u57a3O/6C1h84Kz9vDz3YpJICfL0kSbviL+iRyC36+zyFjsHl9VG/f9h9nbmbjipy3WElXkpTgK+XItrV1FNtajriKwAAUGQw4ganalq1tJpWLZ2rff72+FyhTZJW/3ZGCclXVdHfx+ZrLNgRr5d/iLVuJ15K0yvL4uTrVUyP3V21QHUDAFAUcT8JLpF4KS3PdsOQzl1Kt25vP3JOz3y2U11mbNDYb37RoTOXch0zd9PRPM+VXzsAAGbFiBscLvVappbtSVDMiSRVLu2jXs2qqFxJrxx9WtYsq7X7z+Y6tmwJT9UO9JUkLf/1lIZ8GW0dmfvt1EUt//WUFg8JVe3AktZj/ki6mmcdJ/NpBwDArBhxg0OlpF5Tr8gtGv3NL/p823G9+uNvun/6Ov1yIilHvydaVlWd/wW06ywW6YWu9eRVLPu1WG+u3J/rdurFtAzNWnsoR1vjfN7E0LgKb2gAANxeGHGD3T7fdlxfbD2us5fS9I/qpTWiYx0FV8ye6PHRhiPa90dKjv4XUzM06bt9Wjq0tbXNz9tD3w4O1Vfb47Xl8DkFlPDU4y2r6h/Vy0iSLqVl5HlbVJJ++T05x/aIjnW089h2pf9l7ThPdzeN7FTHId8XAICiguAGu7zz00G9/dMB6/aKfae15dA5/TD8XlUrW0Jr9p/J87iYE0k6fzldaBMAKQAAG/lJREFUZUp4Wtv8vD00oHV1Pdq8isqU8JTl/7d392FRVvn/wN83A4w8DoEogQSiYOKIhH41Uh5sQxK/pluta6uisVkqPqVm6n5dTUolL9pWEx/yq7JmD+5itRlf0h9oGRrKCKWIigJCCiL4wIMPMMP5/cE6OQ7IAKPD4Pt1XXNdzLnPOfdn5rqRj+fc59ySpD1mYyWDs501rtTW6fXl7tRF5/3gns741/RgbP6hAPmXauDb3R6vh/oggCNuRETUyTBxI4PdqFNjy8ECvfLq22psyyjC8hf6wdaq6UvK0kKC3PK3mfl6TQPWfHcan2YWo+a2Gj5d7bDw+T54Xvk4AEBmIWFysLdOknhHzNCeAIBTZVXYcOAcjl+4jiecbfHaMB8M8+1qjI9KRETUIfEeNzLYxWs3UX1b3eSxU2WN06MvDfRo8nhkPzfYyX9L6lam5GHzDwWo+U9/BRW1mLHzGLKKrmjrzHq2N2Y/2xsO8sZ73uzllpgw5AmE+rkir7QKY9dn4Ouciyi4XIsDpy9j0tZMfHakGEUVtUb5vERERB0NR9zIYG4KG9hay3CjTqN3rJdr40KDcYM8cfJiFXb8dF67sGCQ12OIG6vU1r1Rp8bnR0r0+mgQwLZDRRj0n/vcLCwkvBjUA19klaD6tgY1t9XYmVmM7OJruK3W4Fa97vNQhQAW7z4OAHC1l+O/A9zwWmgveDgZviccERFRR8bEjQxmL7fEpGAvbPped7rUxkqGV4d6AwAkScI7Y5R4PawXjv96De5ONnr3ml2prcPNev3kDwAuXNXdwmPFnpO4VKW759vJUt3FD025XHMb2w6dR9Lh81gS1ZdPUSAiok6BU6XUKr8P9IDdf6Yu7+jlaocnnO0AAEIIfH/mMlJ+KYXcSgalu0KvDzfHLuh2z75udwTetbVHnboBB5pZ7GCoBgG8+20ezlyqblc/REREHQFH3KhV/uerE6i9rTtaduJiFV7ZfBijA92x86di5N+1jUdADwX+ETMYTra/rSa1lFnATdEF5dW6I2mSBIwf7AkAqL5Vj29+vti40tQIj9P99pdS+EU4tFyRiIioA2PiRgYrr76FrPNXmzymKr4GVfE1vfJffr2ONd+dxnu/768tK6yoxS/37MUGNOZnu1W/YlSAOyZvO4JrN+qNFrswQvJHRERkapwqJYNZ3LXPWmt8e7xU533+faYtP/6xEFOMnLQB0G4zQkREZM6YuJHButrL4dPVrtXt7h3scn+s+VWeQgBX25G0WUjA3emlJAHzI/zg7+7Y5j6JiIg6Ck6VUqskjBuA3ycealWbwHueJerv5giZBaBpaKZBOwgAn04dgjOXanBbrcFzfbvDx9W+xXZERETmgCNu1CpPPfEY3rtrTzZD/Pxr475rdzQIgS6Wsvu0aDshgOmfqDD2KQ+8HtqLSRsREXUqHHGjVjt0rrJV9a/dqMcPZy7jbHkNci9W4WptHWqb2MT3bhYStBv4tta1m2p8lX0Bk5/xblsHREREHRQTN2qVS1W39BYbGGLGJ8dQ34pMrEE03p82UumGWc/64l+qX/G/PxYa3P5+CyCIiIjMFadKqVXaupFta5K2O4QA/u9EGWysZFj63/5YGNnH4Lb9e+hv/EtERGTumLhRq8gtH+4lIwSQcqJxhG/G8N5454V+eNyx6acu3GEvt8SLQT0eRnhEREQPFRM3apW+jzvCSta2/dzup6eLbbPH6tS/LT+d/Iw3/m9uaLN1LSQgeXowrGS8tImIqPPhXzdqFYcuVnhhgLtR+7SWWWBqaPMPgbe00E0U7eWWcLK1arJuiG9X9HHjnm1ERNQ5MXGjVjtxocpofT37ZDf8a3owHrvrWab3ulSl+0xTS5kFooO99epJEvBaSPMJIBERkbnjqlJqlYvXbuK0EVdspp8qx6WqW3Cxaz5x83Kxxa16Db7OuYC0vHJYWEgI9nHBG2E++OJICa7drIe7UxcsGtkXIb6uRouNiIioo2HiRq1y8fpNo/eZe7H5ETxnO2s85emEYfHpqKip05anniiDs60V6jWN+8FdvHYL7+45CWdbKwxj8kZERJ2UJMS9T5Ls3KqqqqBQKHD9+nU4OvJeqNYQQuC/3tuHihrjPgD+fiwtJHi52OLc5VqD24T7dUVuaTWu36xHqK8r3n6+D3y7OzzAKImIiB4O3uNGBssuufZQkzYAUDeIViVtAHDgTAUuV99GnboB/y/vEv64+SeUV996QBESERE9PEzcyGCFrUygOoortXX4/EiJqcMgIiJqNyZuZLB6TUPLlTooVVHrnq9KRETUETFxI4N5PGZj6hDaLLvkuqlDICIiajcmbmSwYb27mjqENqu+pTZ1CERERO3GxI0MJkkSHOQyU4fRJgJA2XUuUCAiIvPGxI1a5Y2wXqYOoc3M+R49IiIigIkbtdLMZ31hZ21+o259utvD07n5B9kTERGZAyZu1Gr7F4SjT3d7o/d7z7PkjUYmSVjzhwEPpnMiIqKHiE9OoDarrLmNPb+U4tj5K7C1tsSUod6oUzdgxs5jKLnaukdjSRLgai9HefXtliu3gkySsGFiEEb0czNqv0RERKbAxI2MLi3vEv6clKVXbmkhoZ+7I6pu1eNWfQNKm1gsIKFxIcH9yCRAY8BVO6x3Vywa+SSUHgrDAiciIurgOFVKRve7vt3xRqiPztSnh5MNUuaE4OuZw7B/wXD8+Paz6GpvrddWAHjGxwVK96aT6m4Ocmx7dTBsrFq+dMP7uDJpIyKiTsXS1AFQ57Q4qi8mPu2Fw+cq4WxnjbA+rrCS/ZZsVdbeRkVNXZNtr96sx9exQzHxfzNxpPCKtvwxWytsf3Uw/N0doVoagfm7fkZqbhmaGzPu5tjFqJ+JiIjI1Ew+VZqYmIg1a9agtLQU/fr1w4cffoiQkJBm63///feYN28ecnNz4e7ujoULF2LatGkGn49TpcbRe9G3MIctbQf1sMW/Zg43dRhERERGYdKp0i+++AJz587FX/7yF2RnZyMkJAQjR45EcXFxk/ULCwsRFRWFkJAQZGdnY8mSJZg9ezaSk5MfcuSPrn5/+RbeZpK0AUDWrzfgvehbvPN1jqlDISIiajeTjrgNGTIEQUFB2LBhg7asb9++GDt2LFatWqVX/+2338a///1v5OXlacumTZuGn3/+GYcPHzbonBxxa7vqG3Xov2KfqcNos6LVo0wdAhERUbuYbMStrq4OKpUKI0aM0CkfMWIEDh061GSbw4cP69WPjIxEVlYW6uvrm2xz+/ZtVFVV6byobcw5aQOAUX9PN3UIRERE7WKyxK2iogIajQbdu3fXKe/evTvKysqabFNWVtZkfbVajYqKiibbrFq1CgqFQvvy9PQ0zgcgs3PmUuv2liMiIupoTL4diCTpbpcvhNAra6l+U+V3LF68GNevX9e+SkpK2hnxo8vO/J50pWP5qCdNHQIREVG7mCxx69q1K2Qymd7oWnl5ud6o2h1ubm5N1re0tISLi0uTbeRyORwdHXVe1Da575n3PWIThvYydQhERETtYrLEzdraGgMHDsS+fbr3Te3btw/PPPNMk22Cg4P16u/duxeDBg2ClZXVA4uVfnNq+YiWK3VAXJhARESdgUlXlX7xxReYNGkSNm7ciODgYGzevBkff/wxcnNz4eXlhcWLF+PChQv4xz/+AaBxOxClUok33ngDU6dOxeHDhzFt2jR89tlneOmllww6J1eVEhERkbky6ZMT/vjHP6KyshIrVqxAaWkplEolUlJS4OXlBQAoLS3V2dOtZ8+eSElJwZtvvon169fD3d0da9euNThpIyIiIjJnJn9ywsPGETciIiIyVyZfVUpEREREhmHiRkRERGQmmLgRERERmQkmbkRERERmgokbERERkZlg4kZERERkJpi4EREREZkJJm5EREREZoKJGxEREZGZYOJGREREZCaYuBERERGZCSZuRERERGbC0tQBPGxCCACND5snIiIyNgcHB0iSZOowqJN65BK36upqAICnp6eJIyEios7o+vXrcHR0NHUY1ElJ4s4Q1COioaEBFy9e5P+IOpiqqip4enqipKSE/+ARtYC/Lx0b/77Qg/TIjbhZWFigR48epg6DmuHo6Mg/REQG4u8L0aOHixOIiIiIzAQTNyIiIiIzwcSNOgS5XI5ly5ZBLpebOhSiDo+/L0SPrkducQIRERGRueKIGxEREZGZYOJGREREZCaYuFGnFx4ejrlz55o6DKKHoqioCJIkIScnx9ShENEDwMTtETVlyhRIkoTVq1frlH/11VedbuPI3bt3Iy4uztRh0COsrKwMs2bNgo+PD+RyOTw9PTF69GikpaUZ/Vyenp4oLS2FUqk0et9EZHpM3B5hXbp0QXx8PK5evWrqUNqtvr6+2TJnZ2c4ODi0uW+NRoOGhoY2t6dHW1FREQYOHIj09HS8//77OH78OFJTUzF8+HDExsa2qc/mrsm6ujrIZDK4ubnB0rLt+6vX1dW1uS0RPVhM3B5hzz33HNzc3LBq1apm6yQnJ6Nfv36Qy+Xw9vZGQkKCznFvb2+sXLkSMTExcHBwwBNPPIHNmze3eO7c3FyMGjUKjo6OcHBwQEhICM6dOwcAOHr0KCIiItC1a1coFAqEhYXh2LFjOu0lScLGjRsxZswY2NnZ4d1338Xy5csRGBiIrVu3akc2hBB6U6V1dXVYuHAhPDw8YGdnhyFDhuDAgQPa49u3b4eTkxP27NkDf39/yOVynD9/3pCvlEjPjBkzIEkSjhw5gpdffhl+fn7o168f5s2bh59++gkA8MEHH6B///6ws7ODp6cnZsyYgZqaGm0fzV2T3t7eePfddzFlyhQoFApMnTq1yanSkydPIioqCvb29ujevTsmTZqEiooK7fHw8HDMnDkT8+bNQ9euXREREfHwviAiahUmbo8wmUyGlStXYt26dfj111/1jqtUKowbNw7jx4/H8ePHsXz5cixduhTbt2/XqZeQkIBBgwYhOzsbM2bMwPTp03Hq1Klmz3vhwgWEhoaiS5cuSE9Ph0qlQkxMDNRqNQCguroakydPxsGDB/HTTz/B19cXUVFRqK6u1uln2bJlGDNmDI4fP46YmBgAwNmzZ7Fr1y4kJyc3e4/Pq6++ioyMDHz++ef45Zdf8Ic//AHPP/888vPztXVu3LiBVatWYcuWLcjNzUW3bt0M+k6J7nblyhWkpqYiNjYWdnZ2esednJwAND6Kb+3atThx4gSSkpKQnp6OhQsX6tRt7ppcs2YNlEolVCoVli5dqneO0tJShIWFITAwEFlZWUhNTcWlS5cwbtw4nXpJSUmwtLRERkYGNm3aZKyvgIiMTdAjafLkyWLMmDFCCCGefvppERMTI4QQ4ssvvxR3Los//elPIiIiQqfdW2+9Jfz9/bXvvby8xMSJE7XvGxoaRLdu3cSGDRuaPffixYtFz549RV1dnUGxqtVq4eDgIL755httGQAxd+5cnXrLli0TVlZWory8XKc8LCxMzJkzRwghxNmzZ4UkSeLChQs6dX73u9+JxYsXCyGE2LZtmwAgcnJyDIqPqDmZmZkCgNi9e3er2u3atUu4uLho3zd3TXp5eYmxY8fqlBUWFgoAIjs7WwghxNKlS8WIESN06pSUlAgA4vTp00KIxt+RwMDAVsVIRKbBETdCfHw8kpKScPLkSZ3yvLw8DB06VKds6NChyM/Ph0aj0ZYFBARof5YkCW5ubigvLwcAjBw5Evb29rC3t0e/fv0AADk5OQgJCYGVlVWT8ZSXl2PatGnw8/ODQqGAQqFATU0NiouLdeoNGjRIr62XlxdcXV2b/azHjh2DEAJ+fn7auOzt7fH9999rp2oBwNraWudzEbWF+M/+5i0t+Nm/fz8iIiLg4eEBBwcHREdHo7KyErW1tdo6zV2TTf0e3E2lUmH//v061/uTTz4JADrXfEv9EFHH0Pa7V6nTCA0NRWRkJJYsWYIpU6Zoy4UQen9wRBMP2rg3AZMkSXvj9JYtW3Dz5k2dejY2NveNZ8qUKbh8+TI+/PBDeHl5QS6XIzg4WO+G6aamnpoqu1tDQwNkMhlUKhVkMpnOMXt7e+3PNjY2nW51LT18vr6+kCQJeXl5GDt2bJN1zp8/j6ioKEybNg1xcXFwdnbGjz/+iD//+c86i26auyYNueZHjx6N+Ph4vWOPP/64wf0QUcfAxI0AAKtXr0ZgYCD8/Py0Zf7+/vjxxx916h06dAh+fn56SU9zPDw89MoCAgKQlJSE+vr6JkfdDh48iMTERERFRQEASkpKdG6kbo+nnnoKGo0G5eXlCAkJMUqfRM1xdnZGZGQk1q9fj9mzZ+slR9euXUNWVhbUajUSEhJgYdE4CbJr1y6jxRAUFITk5GR4e3u3a6UpEXUMnColAED//v0xYcIErFu3Tls2f/58pKWlIS4uDmfOnEFSUhI++ugjLFiwoF3nmjlzJqqqqjB+/HhkZWUhPz8fO3bswOnTpwEAvXv3xo4dO5CXl4fMzExMmDChxVE6Q/n5+WHChAmIjo7G7t27UVhYiKNHjyI+Ph4pKSlGOQfR3RITE6HRaDB48GAkJycjPz8feXl5WLt2LYKDg9GrVy+o1WqsW7cOBQUF2LFjBzZu3Gi088fGxuLKlSt45ZVXcOTIERQUFGDv3r2IiYnRueWBiMwDEzfSiouL05kKDQoKwq5du/D5559DqVTir3/9K1asWKEzndoWLi4uSE9PR01NDcLCwjBw4EB8/PHH2tG3rVu34urVq3jqqacwadIkzJ4926irOrdt24bo6GjMnz8fffr0wQsvvIDMzEx4enoa7RxEd/Ts2RPHjh3D8OHDMX/+fCiVSkRERCAtLQ0bNmxAYGAgPvjgA8THx0OpVGLnzp333aKntdzd3ZGRkQGNRoPIyEgolUrMmTMHCoVCO8JHROZDEk3dtEREREREHQ7/u0VERERkJpi4EREREZkJJm5EREREZoKJGxEREZGZYOJGREREZCaYuBERERGZCSZuRERERGaCiRsRERGRmWDiRtRJFRUVQZIk5OTkmDoUIiIyEiZuRA9JWVkZZs2aBR8fH8jlcnh6emL06NFIS0t7IOfz9PREaWkplEql0fuura3F22+/DR8fH3Tp0gWurq4IDw/Hnj17tHW8vb3x4Ycftrrv8PBwzJ0715jhEhF1GpamDoDoUVBUVIShQ4fCyckJ77//PgICAlBfX4/vvvsOsbGxOHXqVJv61Wg0kCRJ75mTdXV1sLa2hpubW7vivtPPvaZNm4YjR47go48+gr+/PyorK3Ho0CFUVla263xERNQCQUQP3MiRI4WHh4eoqanRO3b16lXtzwkJCUKpVApbW1vRo0cPMX36dFFdXa09vm3bNqFQKMQ333wj+vbtK2QymSgoKBBeXl4iLi5OTJ48WTg6Ooro6GhRWFgoAIjs7Gxt+9zcXDFy5EhhZ2cnunXrJiZOnCguX76sPR4WFiZiY2PFm2++KVxcXERoaGiTn0ehUIjt27c3+3nDwsIEAJ2XEEJUVFSI8ePHCw8PD2FjYyOUSqX49NNPte0mT56s166wsFD7ue/25Zdfirv/CcvJyRHh4eHC3t5eODg4iKCgIHH06NFmYyQiMkecKiV6wK5cuYLU1FTExsbCzs5O77iTk5P2ZwsLC6xduxYnTpxAUlIS0tPTsXDhQp36N27cwKpVq7Blyxbk5uaiW7duAIA1a9ZAqVRCpVJh6dKleucpLS1FWFgYAgMDkZWVhdTUVFy6dAnjxo3TqZeUlARLS0tkZGRg06ZNTX4mNzc3pKSkoLq6usnju3fvRo8ePbBixQqUlpaitLQUAHDr1i0MHDgQe/bswYkTJ/D6669j0qRJyMzMBAD8/e9/R3BwMKZOnapt5+np2dxXq2PChAno0aMHjh49CpVKhUWLFsHKysqgtkRE5oJTpUQP2NmzZyGEwJNPPtli3bvv7erZsyfi4uIwffp0JCYmasvr6+uRmJiIAQMG6LR99tlnsWDBAu37oqIineMbNmxAUFAQVq5cqS3bunUrPD09cebMGfj5+QEAevfujffff/++cW7evBkTJkyAi4sLBgwYgGHDhuHll1/G0KFDAQDOzs6QyWRwcHDQma718PDQiXHWrFlITU3FP//5TwwZMgQKhQLW1tawtbVt9TRvcXEx3nrrLe337Ovr26r2RETmgCNuRA+YEAIAIElSi3X379+PiIgIeHh4wMHBAdHR0aisrERtba22jrW1NQICAvTaDho06L59q1Qq7N+/H/b29trXnSTn3LlzBvcDAKGhoSgoKEBaWhpeeukl5ObmIiQkBHFxcfdtp9Fo8N577yEgIAAuLi6wt7fH3r17UVxc3OI5WzJv3jy89tpreO6557B69Wqdz0RE1FkwcSN6wHx9fSFJEvLy8u5b7/z584iKioJSqURycjJUKhXWr18PoHGU7Q4bG5smk8CmpmHv1tDQgNGjRyMnJ0fnlZ+fj9DQUIP7ucPKygohISFYtGgR9u7dixUrViAuLg51dXXNtklISMDf/vY3LFy4EOnp6cjJyUFkZOR92wCNU8h3EuA77v5OAGD58uXIzc3FqFGjkJ6eDn9/f3z55ZcGfRYiInPBxI3oAXN2dkZkZCTWr1+vM3J2x7Vr1wAAWVlZUKvVSEhIwNNPPw0/Pz9cvHjRaHEEBQUhNzcX3t7e6N27t87L0GTtfvz9/aFWq3Hr1i0AjSODGo1Gp87BgwcxZswYTJw4EQMGDICPjw/y8/N16jTVztXVFdXV1TrfX1P70/n5+eHNN9/E3r178eKLL2Lbtm3t/lxERB0JEzeihyAxMREajQaDBw9GcnIy8vPzkZeXh7Vr1yI4OBgA0KtXL6jVaqxbtw4FBQXYsWMHNm7caLQYYmNjceXKFbzyyis4cuQICgoKsHfvXsTExOglSi0JDw/Hpk2boFKpUFRUhJSUFCxZsgTDhw+Ho6MjgMZ93H744QdcuHABFRUVABrvn9u3bx8OHTqEvLw8vPHGGygrK9Pp29vbG5mZmSgqKkJFRQUaGhowZMgQ2NraYsmSJTh79iw+/fRTbN++Xdvm5s2bmDlzJg4cOIDz588jIyMDR48eRd++fdv3pRERdTBM3Igegp49e+LYsWMYPnw45s+fD6VSiYiICKSlpWHDhg0AgMDAQHzwwQeIj4+HUqnEzp07sWrVKqPF4O7ujoyMDGg0GkRGRkKpVGLOnDlQKBR6+8C1JDIyEklJSRgxYgT69u2LWbNmITIyErt27dLWWbFiBYqKitCrVy+4uroCAJYuXYqgoCBERkYiPDwcbm5uGDt2rE7fCxYsgEwmg7+/P1xdXVFcXAxnZ2d88sknSElJQf/+/fHZZ59h+fLl2jYymQyVlZWIjo6Gn58fxo0bh5EjR+Kdd95p+xdGRNQBSeLeG0eIiIiIqEPiiBsRERGRmWDiRkRERGQmmLgRERERmQkmbkRERERmgokbERERkZlg4kZERERkJpi4EREREZkJJm5EREREZoKJGxEREZGZYOJGREREZCaYuBERERGZCSZuRERERGbi/wOc2vQZua7+oAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 644x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = sns.catplot(data=misexp_exp_sv_carriers_df, \n",
    "           x=\"status\", \n",
    "            y=\"TPM\", \n",
    "                 hue=\"misexp\"\n",
    "           )\n",
    "ax.set(ylabel=\"ROPN1B Expression (TPM)\", xlabel=\"Carrier Status\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16e39136",
   "metadata": {},
   "source": [
    "### Transcript quantification data (Salmon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "dae4a686",
   "metadata": {},
   "outputs": [],
   "source": [
    "# misexpressed genes transcripts \n",
    "gencode_path = wkdir_path.joinpath(\"reference/gencode/gencode.v31.annotation.sorted.gtf.gz\")\n",
    "misexp_gene_transcript_ids = []\n",
    "for gtf in pysam.TabixFile(str(gencode_path)).fetch(parser = pysam.asGTF()): \n",
    "    if gtf.gene_id.split(\".\")[0] == misexp_gene_id and gtf.feature == \"transcript\": \n",
    "        misexp_gene_transcript_ids.append(gtf.transcript_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "328a1ddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# collect transcript quantification files\n",
    "quant_sf_dir = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq_n5188_v97/results/salmon\"\n",
    "quant_sf_dir_path = Path(quant_sf_dir)\n",
    "misexp_transcript_df = pd.DataFrame(columns=[\"Name\"])\n",
    "for rna_id in rna_id_carriers: \n",
    "    rna_id_path = f\"{rna_id}.quant.sf\"\n",
    "    quant_sf_path = quant_sf_dir_path.joinpath(rna_id_path)\n",
    "    quant_sf_df = pd.read_csv(quant_sf_path, sep=\"\\t\")\n",
    "    quant_sf_df = quant_sf_df.rename(columns={\"TPM\": rna_id})\n",
    "    misexp_quant_sf_df = quant_sf_df[quant_sf_df.Name.isin(misexp_gene_transcript_ids)]\n",
    "    misexp_transcript_df = pd.merge(misexp_transcript_df, \n",
    "                                    misexp_quant_sf_df[[\"Name\", rna_id]], \n",
    "                                    on=\"Name\", \n",
    "                                    how=\"outer\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "717460d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "misexp_transcript_tp_df = misexp_transcript_df.set_index(\"Name\").transpose()\n",
    "misexp_transcript_tp_proport_df = misexp_transcript_tp_df.div(misexp_transcript_tp_df.sum(axis=1), axis=0) * 100\n",
    "misexp_transcript_tp_proport_df = misexp_transcript_tp_proport_df.reset_index().melt(id_vars=\"index\").rename(columns={\"index\": \"rna_id\", \"Name\": \"transcript_id\", \"value\": \"perc_expression\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "dff09c8b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of RNA IDs: 10\n"
     ]
    }
   ],
   "source": [
    "print(f\"Number of RNA IDs: {len(misexp_transcript_tp_proport_df.rna_id.unique())}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "2c3a715c",
   "metadata": {},
   "outputs": [],
   "source": [
    "misexp_transcript_tp_proport_path = out_dir.joinpath(\"ropn1b_proport_expression.tsv\")\n",
    "misexp_transcript_tp_proport_df.to_csv(misexp_transcript_tp_proport_path, sep=\"\\t\", index=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
